clear
rng default



%% Empirical distribution from data on cohort 1950-1951

%%Men%%
year=1950;
% Probability distribution of Ychosen|X
condprobsmen_temp=load('probs_allyears_data.mat',['condprobsmen' num2str(year) '_140']);
condprobsmen=condprobsmen_temp.(['condprobsmen' num2str(year) '_140']);
PYX_cond_emp=[condprobsmen(:,2:end) condprobsmen(:,1)]; %5x6
%Order: P(Ychosen=1|X=1) P(Ychosen=2|X=1) P(Ychosen=3|X=1) P(Ychosen=4|X=1) P(Ychosen=5|X=1) P(Ychosen=0|X=1) 
%       P(Ychosen=1|X=2) P(Ychosen=2|X=2) P(Ychosen=3|X=2) P(Ychosen=4|X=2) P(Ychosen=5|X=2) P(Ychosen=0|X=2) 
%       P(Ychosen=1|X=3) P(Ychosen=2|X=3) P(Ychosen=3|X=3) P(Ychosen=4|X=3) P(Ychosen=5|X=3) P(Ychosen=0|X=3) 
%       P(Ychosen=1|X=4) P(Ychosen=2|X=4) P(Ychosen=3|X=4) P(Ychosen=4|X=4) P(Ychosen=5|X=4) P(Ychosen=0|X=4) 
%       P(Ychosen=1|X=5) P(Ychosen=2|X=5) P(Ychosen=3|X=5) P(Ychosen=4|X=5) P(Ychosen=5|X=5) P(Ychosen=0|X=5) 
probsmen_temp=load('probs_allyears_data.mat',['probsmen' num2str(year) '_140']);
PX_emp=probsmen_temp.(['probsmen' num2str(year) '_140']); %5x1
%Order: P(X=1) P(X=2) P(X=3) P(X=4) P(X=5) 

%%Women%%
year=1951;
% Probability distribution of Xchosen|Y
condprobswomen_temp=load('probs_allyears_data.mat',['condprobswomen' num2str(year) '_140']);
condprobswomen=condprobswomen_temp.(['condprobswomen' num2str(year) '_140']);
PXY_cond_emp=[condprobswomen(:,2:end) condprobswomen(:,1)]; %5x6
%Order: P(Xchosen=1|Y=1) P(Xchosen=2|Y=1) P(Xchosen=3|Y=1) P(Xchosen=4|Y=1) P(Xchosen=5|Y=1) P(Xchosen=0|Y=1) 
%       P(Xchosen=1|Y=2) P(Xchosen=2|Y=2) P(Xchosen=3|Y=2) P(Xchosen=4|Y=2) P(Xchosen=5|Y=2) P(Xchosen=0|Y=2) 
%       P(Xchosen=1|Y=3) P(Xchosen=2|Y=3) P(Xchosen=3|Y=3) P(Xchosen=4|Y=3) P(Xchosen=5|Y=3) P(Xchosen=0|Y=3) 
%       P(Xchosen=1|Y=4) P(Xchosen=2|Y=4) P(Xchosen=3|Y=4) P(Xchosen=4|Y=4) P(Xchosen=5|Y=4) P(Xchosen=0|Y=4) 
%       P(Xchosen=1|Y=5) P(Xchosen=2|Y=5) P(Xchosen=3|Y=5) P(Xchosen=4|Y=5) P(Xchosen=5|Y=5) P(Xchosen=0|Y=5) 
probswomen_temp=load('probs_allyears_data.mat',['probswomen' num2str(year) '_140']);
PY_emp=probswomen_temp.(['probswomen' num2str(year) '_140']); %5x1
%Order: P(Y=1) P(Y=2) P(Y=3) P(Y=4) P(Y=5) 


%% Get the systematic match surplus 
%Notation: u12=man is type 2 and chooses woman of type 1
ntypes_m=5;
ntypes_w=5;
u1_emp=zeros(1,ntypes_m);
u2_emp=zeros(1,ntypes_m);
u3_emp=zeros(1,ntypes_m);
u4_emp=zeros(1,ntypes_m);
u5_emp=zeros(1,ntypes_m);
for x=1:ntypes_m
    u1_emp(x)=log(PYX_cond_emp(x,1)/PYX_cond_emp(x,6)); %u1(x)=log(P(Ychosen=1|X=x)/P(Ychosen=0|X=x))
    u2_emp(x)=log(PYX_cond_emp(x,2)/PYX_cond_emp(x,6)); %u2(x)=log(P(Ychosen=2|X=x)/P(Ychosen=0|X=x))
    u3_emp(x)=log(PYX_cond_emp(x,3)/PYX_cond_emp(x,6)); %u3(x)=log(P(Ychosen=3|X=x)/P(Ychosen=0|X=x))
    u4_emp(x)=log(PYX_cond_emp(x,4)/PYX_cond_emp(x,6)); %u4(x)=log(P(Ychosen=4|X=x)/P(Ychosen=0|X=x))
    u5_emp(x)=log(PYX_cond_emp(x,5)/PYX_cond_emp(x,6)); %u5(x)=log(P(Ychosen=5|X=x)/P(Ychosen=0|X=x))
end

v1_emp=zeros(1,ntypes_w);
v2_emp=zeros(1,ntypes_w);
v3_emp=zeros(1,ntypes_w);
v4_emp=zeros(1,ntypes_w);
v5_emp=zeros(1,ntypes_w);
for y=1:ntypes_w
    v1_emp(y)=log(PXY_cond_emp(y,1)/PXY_cond_emp(y,6)); %v1(x)=log(P(Xchosen=1|Y=y)/P(Xchosen=0|Y=y))
    v2_emp(y)=log(PXY_cond_emp(y,2)/PXY_cond_emp(y,6)); %v2(x)=log(P(Xchosen=2|Y=y)/P(Xchosen=0|Y=y))
    v3_emp(y)=log(PXY_cond_emp(y,3)/PXY_cond_emp(y,6)); %v3(x)=log(P(Xchosen=2|Y=y)/P(Xchosen=0|Y=y))
    v4_emp(y)=log(PXY_cond_emp(y,4)/PXY_cond_emp(y,6)); %v4(x)=log(P(Xchosen=2|Y=y)/P(Xchosen=0|Y=y))
    v5_emp(y)=log(PXY_cond_emp(y,5)/PXY_cond_emp(y,6)); %v5(x)=log(P(Xchosen=2|Y=y)/P(Xchosen=0|Y=y))
end


u_emp=[u1_emp(1) u2_emp(1) u3_emp(1) u4_emp(1) u5_emp(1); ...
       u1_emp(2) u2_emp(2) u3_emp(2) u4_emp(2) u5_emp(2); ...
       u1_emp(3) u2_emp(3) u3_emp(3) u4_emp(3) u5_emp(3); ...
       u1_emp(4) u2_emp(4) u3_emp(4) u4_emp(4) u5_emp(4); ...
       u1_emp(5) u2_emp(5) u3_emp(5) u4_emp(5) u5_emp(5)];

v_emp=[v1_emp(1) v1_emp(2) v1_emp(3) v1_emp(4) v1_emp(5); ...
       v2_emp(1) v2_emp(2) v2_emp(3) v2_emp(4) v2_emp(5); ...
       v3_emp(1) v3_emp(2) v3_emp(3) v3_emp(4) v3_emp(5); ...
       v4_emp(1) v4_emp(2) v4_emp(3) v4_emp(4) v4_emp(5); ...
       v5_emp(1) v5_emp(2) v5_emp(3) v5_emp(4) v5_emp(5)];

Phi_emp=u_emp+v_emp;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Primitives
nm=6000; %number of men
nw=6000; %number of women
supX=[1 2 3 4 5]; %support of X (men's types)
supY=[1 2 3 4 5]; %support of Y (women's types)
supXY=sup(supX, supY); %support of (X,Y)
%supXY order:
%(1,1)
% ...
%(5,1)
%(1,2)
% ...
%(5,2)
%(1,3)
% ...
%(5,3)
%(1,4)
% ...
%(5,4)
%(1,5)
% ...
%(5,5)


PhiXY=Phi_emp;
%PhiXY=
%Phi11 Phi12 ... Phi15
%Phi21 Phi22 ... Phi25
%...
%Phi51 Phi52 ... Phi55

%% Generate data
%Men
X=randsmpl(PX_emp, nm,1); %nmx1;
epsilon = -evrnd(0,1,nm,ntypes_w+1); %nmxntypes_w+1; 
%epsilon=(epsilon1, epsilon2, ..., epsilon0) i.i.d. Gumbel with location 0 and scale 1, independent of X
%Gumbel defined as in Wikipedia (Matlab flips the signs)

%Women
Y=randsmpl(PY_emp, nw,1); %nwx1;
eta = -evrnd(0,1,nw,ntypes_m+1); %nwxntypes_m+1; 
%eta=(eta1, eta2, ..., eta0) i.i.d. Gumbel with location 0 and scale 1, independent of Y

%% Equilibrium
%(1) Compute total surplus for each possible couple
Phi=zeros(nm,nw); %nmxnw
for j=1:nw %for each woman
    for i=1:nm %for each man
        Phi(i,j)=PhiXY(X(i), Y(j))+epsilon(i,Y(j))+eta(j,X(i));
    end
end

%(2) Reshape Phi into a vector (nm*nw)x1 following the order
%man1 woman1
%man2 woman1
%...
%man1 woman2
%man2 woman2
%...
%man1 woman3
%man2 woman3
%...
Phir=reshape(Phi,nm*nw,1); %(nm*nw)x1

%(3) Augment Phir for singles following the order
%man1 woman1
%man1 woman2
%man1 woman3
%...
%man2 woman1
%man2 woman2
%man2 woman3
%...
%man1 single
%man2 single
%...
%woman1 single
%woman2 single
%woman3 single
%...

Phira=[Phir; epsilon(:,end); eta(:,end)]; %(nm*nw+nm+nw)x1

%(3) Find equilibrium matching 
mu=equilibrium(nw, nm,Phira); %(nm*nw+nm+nw)x1


%% Conditional probabilities

% Men
Ychosen=zeros(nm,1); %nmx1 listing the type of woman man i is matched with
for i=1:nm
    %select from mu the (nw+1)x1 sub-vector reporting the choice of man i
    index=(1:nm:(nm*nw+1)).'+(i-1);
    submu=mu(index);
    %save the type of the woman man i is matched with
    if submu(end)==0 %man i has not chosen to be single
        Ychosen(i)=Y(logical(submu(1:end-1)));
    else %man i has chosen to be single
        Ychosen(i)=0;
    end
end

% Probability distribution of Ychosen|X
PX=zeros(ntypes_m,1);
for x=1: ntypes_m
    PX(x)=sum(X==x)/nm;
end
PYX=zeros(ntypes_m, ntypes_w+1);
for x=1:ntypes_m
    for y=1:ntypes_w
        PYX(x,y)=sum(X==x & Ychosen==y)/nm;
    end
end
for x=1:ntypes_m
    PYX(x,ntypes_w+1)=sum(X==x & Ychosen==0)/nm;
end 
PYX_cond=PYX./kron(PX,ones(1,ntypes_w+1)); 
%order: P(Ychosen=1|X=1) P(Ychosen=2|X=1) P(Ychosen=3|X=1) P(Ychosen=4|X=1) P(Ychosen=5|X=1) P(Ychosen=0|X=1) 
%       ...
%       P(Ychosen=1|X=5) P(Ychosen=2|X=5) P(Ychosen=3|X=5) P(Ychosen=4|X=5) P(Ychosen=5|X=5) P(Ychosen=0|X=5) 


% Women
Xchosen=zeros(nw,1); %nwx1 listing the type of man woman j is matched with
for j=1:nw
    %select from mu the (nm+1)x1 sub-vector reporting the choice of woman j
    index=[(1:nm).'+(j-1)*nm; nm*nw+nm + j];
    submu=mu(index);
    %save the type of the man woman j is matched with
    if submu(end)==0 %woman j has not chosen to be single
        Xchosen(j)=X(logical(submu(1:end-1)));
    else %woman j has chosen to be single
        Xchosen(j)=0;
    end
end

% Probability distribution of Xchosen|Y
PY=zeros(ntypes_w,1);
for y=1: ntypes_w
    PY(y)=sum(Y==y)/nw;
end
PXY=zeros(ntypes_w, ntypes_m+1);
for y=1:ntypes_w
    for x=1:ntypes_m
        PXY(y,x)=sum(Y==y & Xchosen==x)/nw;
    end
end
for y=1:ntypes_w
    PXY(y,ntypes_m+1)=sum(Y==y & Xchosen==0)/nw;
end 
PXY_cond=PXY./kron(PY,ones(1,ntypes_m+1)); 
%order: P(Xchosen=1|Y=1) P(Xchosen=2|Y=1) P(Xchosen=3|Y=1) P(Xchosen=4|Y=1) P(Xchosen=5|Y=1) P(Xchosen=0|Y=1)
%       ...
%       P(Xchosen=1|Y=5) P(Xchosen=2|Y=5) P(Xchosen=3|Y=5) P(Xchosen=4|Y=5) P(Xchosen=5|Y=5) P(Xchosen=0|Y=5) 


save('PYX_cond.mat', 'PYX_cond')
save('PXY_cond.mat', 'PXY_cond')

%% Systematic payoff  
u1=zeros(1,ntypes_m);
u2=zeros(1,ntypes_m);
u3=zeros(1,ntypes_m);
u4=zeros(1,ntypes_m);
u5=zeros(1,ntypes_m);
for x=1:ntypes_m
    u1(x)=log(PYX_cond(x,1)/PYX_cond(x,6)); %u1(x)=log(P(Ychosen=1|X=x)/P(Ychosen=0|X=x))
    u2(x)=log(PYX_cond(x,2)/PYX_cond(x,6)); %u2(x)=log(P(Ychosen=2|X=x)/P(Ychosen=0|X=x))
    u3(x)=log(PYX_cond(x,3)/PYX_cond(x,6)); %u3(x)=log(P(Ychosen=3|X=x)/P(Ychosen=0|X=x))
    u4(x)=log(PYX_cond(x,4)/PYX_cond(x,6)); %u4(x)=log(P(Ychosen=4|X=x)/P(Ychosen=0|X=x))
    u5(x)=log(PYX_cond(x,5)/PYX_cond(x,6)); %u5(x)=log(P(Ychosen=5|X=x)/P(Ychosen=0|X=x))
end

v1=zeros(1,ntypes_w);
v2=zeros(1,ntypes_w);
v3=zeros(1,ntypes_w);
v4=zeros(1,ntypes_w);
v5=zeros(1,ntypes_w);
for y=1:ntypes_w
    v1(y)=log(PXY_cond(y,1)/PXY_cond(y,6)); %v1(x)=log(P(Xchosen=1|Y=y)/P(Xchosen=0|Y=y))
    v2(y)=log(PXY_cond(y,2)/PXY_cond(y,6)); %v2(x)=log(P(Xchosen=2|Y=y)/P(Xchosen=0|Y=y))
    v3(y)=log(PXY_cond(y,3)/PXY_cond(y,6)); %v3(x)=log(P(Xchosen=2|Y=y)/P(Xchosen=0|Y=y))
    v4(y)=log(PXY_cond(y,4)/PXY_cond(y,6)); %v4(x)=log(P(Xchosen=2|Y=y)/P(Xchosen=0|Y=y))
    v5(y)=log(PXY_cond(y,5)/PXY_cond(y,6)); %v5(x)=log(P(Xchosen=2|Y=y)/P(Xchosen=0|Y=y))
end


u=[u1(1) u2(1) u3(1) u4(1) u5(1); ...
   u1(2) u2(2) u3(2) u4(2) u5(2); ...
   u1(3) u2(3) u3(3) u4(3) u5(3); ...
   u1(4) u2(4) u3(4) u4(4) u5(4); ...
   u1(5) u2(5) u3(5) u4(5) u5(5)];

v=[v1(1) v1(2) v1(3) v1(4) v1(5); ...
   v2(1) v2(2) v2(3) v2(4) v2(5); ...
   v3(1) v3(2) v3(3) v3(4) v3(5); ...
   v4(1) v4(2) v4(3) v4(4) v4(5); ...
   v5(1) v5(2) v5(3) v5(4) v5(5)];


U_CS=zeros(ntypes_m,ntypes_w);
for i=1:ntypes_m
for j=1:ntypes_w
    U_CS(i,j)=u(i,j);  
end
end

V_CS=zeros(ntypes_m,ntypes_w);
for j=1:ntypes_w
for i=1:ntypes_m
    V_CS(i,j)=v(i,j);  
end
end


save('U_CS.mat', 'U_CS')
save('V_CS.mat', 'V_CS')




